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The Yukawa model in the quenched approximation is expressed as a disordered statistical me- 
chanics model on a 4-dimensional Euclidean lattice. We study this model. A particular attention 
is given to the singularities of the Dirac operator in the phase diagram. A careful analysis of a 
particular limiting case shows that finite volume effects can be huge and questions the quenched 
approximation. This is confirmed by a Monte-Carlo simulation in this limiting case and without the 
quenched approximation. We include also some results concerning the symmetries of this model. 



Since long it has been recognized that a Quantum Field Theory can be expressed as a Statistical Mechanics problem. 
The two areas have benefited from this proximity and many techniques developed in one context have then been used in 
the other [1] . The expression of quantum field theories as statistical mechanical problems has been especially useful in 
the context of the fundamental theory describing the interaction of quarks and gluons i.e. Quantum ChromoDynamics 
(QCD), where usual perturbative techniques fail. Indeed it has been found that perturbative series in any quantum 
field theory have a zero convergence radius and are asymptotic but never convergent In such situations, it is 
common to resort to a numerical approach based on the Feynman path integral formulation, where the system is 
described by a discretized action on a space-time lattice [3]. 

The numerical formulation of QCD on a lattice is nowadays among the most challenging problems of numerical 
physics and the progress have been very important during the last decades. The methods developed in this context 
can also be applied to other quantum field theories |4] in situations where perturbation theory fails, for example in 
the investigation of binding energies. Indeed, such calculations would require the evaluation of an infinite number of 
contributions in a perturbative scheme. 

In this paper we study the simplest fermion quantum field theory in four space-time dimensions, that is the model 
introduced by H. Yukawa for the nuclear interaction 5 . The model is described in detail in reference ^6,, but since 
the aim of this paper is to adopt a statistical mechanics point of view, we simply sketch extremely schematically how 
one goes from the nuclear physics modelization to the statistical mechanics formulation. 

Similar models were analyzed some time ago using the same techniques 0, and two distinct regimes were found: 
for small and large values of the coupling constant the system was numerically solvable while for intermediate values 
it was not. In this paper we shall address in detail this issue, also found in [B]. The same techniques have been used 
for a numerical study of a similar model W, where some bound on the Higgs boson mass is established based on a 
Yukawa coupling between quarks and the Higgs boson. 



The model introduced by Yukawa aimed at a description of nuclei via the exchange of massive particles in analogy 
with Quantum Electrodynamics, except that the particles mediating the nuclear force have to be massive in order to 
have a finite range interaction. Although we assume QCD is the fundamental theory of quark interactions responsible 
for nuclear interactions, one boson exchange models are still mandatory in the nuclear physics community. 



I. INTRODUCTION AND MODEL 
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Introduction 
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Therefore, in the Yukawa model, nucleons and mesons are considered elementary particles -i.e. without an internal 
structure-, represented by local fields. The mesons are bosons represented by a complex scalar field (f) while nucleons 
are fermions represented by a four component grassmannian Dirac spinor ip. To get a statistical mechanics model one 
works in an Euclidean space instead of a Minkowski space, this is achieved by performing a Wick rotation ^ and the 
space-time is discretized into a four-dimensional hyper-cubic lattice. One possible choice for the discretized action [3] 
is 



S 



+ '^tpxDy,tpx + g^1px(l>X-4'x (1) 



which is the sum of three terms S — Skg + "^w + 5*1. In the first term, which is just a Klein-Gordon action for a free 
bosonic field, x runs over the N sites, v runs over the four space-time direction, and fi is the meson mass. The second 
term is a bilinear in the Dirac- Wilson operator D^. It is a 4iV x AN matrix, with elements 

iDw)xy = ^4Sx,v - ((I4 +liy)Sx,y-iy + (I4 " liy)Sx,y+,y) (2) 

V 

where I4 is a 4 x 4 unity matrix and 7^ are the Dirac matrices ("i/" is the conjugate of '0), k is the so-called hopping 
parameter related to the bare fermion mass M — 1/k — 8. The coupling between the two fields is realized in the 
simplest way by the third term where g is the coupling constant. Every dimensional quantity has been redefined in 
terms of the lattice spacing a, therefore the model depends on the three adimensionalized lattice parameters g, /i and 
K. It depends also on the size of the lattice. In this work we use periodic boundary conditions and take the four 
dimension equal. 

Propagators in Quantum Field Theory are expressed using Wick contractions. From the statistical mechanics point 
of view it amounts to computing expectation values and to combining them together. For example the elementary 
fermion propagator reads: 

S(x, y)^\\ m {mr^)xy det (3) 



where x and y are two sites of the lattice and 



D(0) ^B^ + g^ (4) 



is the interacting Dirac operator. Z is the normalization factor of the probability distribution of the fields and it is 
not calculated in practice. Propagators like ([s]) provide a simple way of computing the renormalized mass m of an 
interacting particle in a QFT, as 

C(x/i) = ^ S'(x,0) '-^ cosh™ — 0:4 j (5) 



X\ ,X2 ,3:3 



where X4 is the time coordinate. The calculation of renormalized masses is performed by producing the fields (px 
according to a joint probability distribution: 

n({0J)~det(7^(0))e-^-«(^) (6) 

and computing S{x,0) as the average over field configurations of (-0(0)^^)^^. Note that it implies solving a linear 
system, not a full inversion of the Dirac operator. 

In the bosons probability distribution Eq. |6] the evaluation of the fermionic determinant is -by large- the most 
expensive part of the calculation. Sophisticated methods have been developed for dealing with this difficulty, as 
Hybrid Monte Carlo simulations ^U\, but the study of the model neglecting the effect of the determinant on the 
weight of field configuration, called quenched approximation, deserves interest yet and will be described in some 
detail in the next section. Finally, let us remind that to extract physical quantities one needs to be as close as possible 
to a critical point so that the operator D{(f)) has low modes. This implies numerical difficulties as in the vicinity of 
any critical point. 



II. THE QUENCHED APPROXIMATION 



The quenched approximation consists in neglecting the variation of det {D{(f))) among the field configurations. From 
a physical point of view, this determinant accounts for the creation of virtual nucleon-anti-nucleon pairs, and its effect 
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is expected to be small as long as meson mass is smaller than nucleon one. It simplifies considerably the problem 
since now Eq. [6] becomes 



Prob({(/)J)^exp(-i^ 



(8 + /i') (bl~2Y,^. 



(7) 



This distribution does not anymore involve the Dirac operator and is easy to implement. Indeed the quadratic form 
in the exponential can be diagonalized straightforwardly, simply going to the discrete Fourier space. We note 4>k the 
Fourier transform of the (j)^- The (j)k are complex and their joint probability factorizes 

prob({,^,.}) ^n^^p l^-^^j 

with cr^ = 2_^_'^ fc" ^ where ki, = 2 sin ^ with the extra constraint — (j)^k in order to get real values for ipx- It 

is then simple to draw independently the real and imaginary part of each (for fc > 0) from a centered Gaussian 
distribution with variance af.. The partial distribution of the {i.e. integrating out all 0y but 0^;) is also a Gaussian 
with a variance a independent of x and given by 

a2(^) = lya2 = ly ^ ^ (9) 

In summary the (j)x are Gaussian dependent with the same variance and the (pk are independent with a variance 
depending on k. 

It is straightforward to compute analytically the meson correlator 

C(t) = ( (t){xi,X2,x:i,Xi)(j){xi + X, X2 + y, X3 + z,X4 + t)\ 

\x,y,z I 

(C does not depend upon the x^''s due to translational invariance) however when the bosons do interact, directly or 
through the fermions when quenched approximation is not assumed, the analytical calculation is not possible and 
one has to perform numerical calculation sampling the field configurations in order to get the re-normalized meson 
mass. In this context, three estimators of the correlator C{t) are possible. For the first estimator Co(i) the point 
(xo, 2/0: -Zo,, io) is fixed and can be chosen to be the origin, for the second estimator Ci{t) the point (xq, J/Oj ^o,j ^o) 
runs over all the points of the time-slice Iq = 0, and for the third one C2{t) all pairs of lattice points are considered. 
Obviously these three estimators give the same average as it should due to the translational invariance, but their 
variances are very different. In appendix [B] we give the three expressions of the variance corresponding to the three 
estimators. We see that only C2(t) is self-averaging. For the first estimator the variance diverges with the size of the 
lattice, while for the second it goes to a finite value. It means that only with the third estimator larger sizes imply 
less configurations in the average. Consequently for the case of interacting bosons, in the unquenched calculation for 
example, the third estimator C2 should be considered. 

III. SYMMETRIES 

In this section we discuss some symmetries of the Dirac operator with Yukawa coupling Eq|4j Being associated to 
the action, these symmetries hold in both quenched and unquenched calculations. They are interesting per se but 
also useful for numerical treatment. 

A. Symmetries holding separately on each boson configuration 

Let us first note that, using the representation for the Dirac matrices the operator J = 17173 is an involution 
verifying J^uJ = transpose (71^) for the four Dirac matrices 71,. It is then straightforward to verify that 

D = JD*J (10) 
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where D* is the complex conjugate of D. Let V be an eigenvector of D belonging to the eigenvalue A. Introducing 
the complex conjugation operator the vector W = JK{V) is also an eigenvector of D belonging to the eigenvalue 
A*. Indeed DJK{V) = JD*K{V) = JK{DV) = JK{\V) = \*JK{V). Moreover V and W are orthogonal. So the 
eigenvalues appear in pair of conjugate values and therefore the determinant is never negative. This non-negativity 
property is useful to perform hybrid Monte-Carlo simulation in the unquenched calculation. 

The relation Eq 10 has another useful consequence. In Eq. [3]S'(x,?/) is a 4 x 4 matrix where row and column are 



indexed by the spin at sites x and y. To compute S{x, y) one solves for the propagator X„ the four linear system 
with the four right hand-side (source term) Y„ 

DX, = (11) 

corresponding to the 4 spin states cr = 0, 1, 2, 3. The 4x4 matrix S{x, y) is obtained selecting the proper elements of 
the four vectors X„. We will show a relation between Xi and X2, obviously this relation hods also between X3 and 



X4. Indeed it is readily verified, using sing Eq. 10 , that D {iJ K{Xi)) — Y2, in other words 

Xi = -iJK{X2) (12) 
and consequently each correlation matrix, and for any field configuration, has the following form 

(a b c d \ 
-/* e* -h* g* J 

and the trace of any of these matrix is simply 2 {Ti{a) + TZ{g)). Note that this form of the correlation matrix holds 
also for any composite particle correlator (even using the so-called smeared source). 

B. Symmetries holding on the average 

We now show that another simplification appears when averaging the correlation matrices over the fields config- 
urations. Let us introduce the automorphy group of the lattice, i.e. permutations tt of the sites of the lattice such 
that the images of two neighboring sites are also two neighboring sites. For any such permutation the two fields 
configurations (px and ipT^i^x) have the same probability, since both the fermions and the bosons actions are invariant 
under the permutation tt. Note that this equality also holds without the quenched approximation. Let's denote tti 
the particular permutation defined by 

'Ki{xi,X2,Xz,Xi) = (-a::i,X2,X3,X4) (14) 

it is clear that tti belongs to the automorphy group. We also introduce 1x2, tts and 7r4 corresponding respectively to 
X2„ X3 and X4. We have 

S{x, (j>x) = 757/c5' (7rfe(x), 7^75 (15) 

where x = (xi, a;2, X3, 0:4) and k — 1,2,3,4, 75 — 71727374 (see ref [Tj section 8.2 for the free fermion case, the 
extension to the Yukawa model treated in this paper is straightforward) . Using this relation the 1-fermion correlation 
matrix, when the source is located at the origin, takes the forms 

c{xi) 

CM= E SM=| » l I (16) 



2:1,2:2,3:3 



c(L4-a;4) 
GO c(L4-a;4) 



the precise form [16] obviously depends on the chosen representation for the Dirac matrices, but in any representation 
the matrix C{t) depends on a single function c{t) instead of 16 functions. 

IV. THE DIRAC OPERATOR SPECTRUM IN THE PHASE SPACE K-g 

Let us recall that the model depends on three independent parameters, k, g and fi. As shown above, in the quenched 
approximation the probability of a (j^x depends only on /i and not on k or g, it is the same everywhere in the parameter 
space. In this section we work at constant value of /i 0.1 
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Any numerical computation of a physical quantity will imply some inversions of the Dirac operator Eq. [4j We know 
that this inversion will have to be performed with values of g and k such that the linear system is difficult to invert. 
In practice, in some region of the <? — k plane and for a given value of the linear sizes of the lattice, solving for X the 
system DX — Y will not be possible. Indeed, depending on the numerical method used, either the algorithm will not 
converge, or it will find a wrong solution. To quantify how ill conditioned the linear system is, it is customary to use 
the condition number. By definition a condition number measure how the solution of the system changes when the 
RHS term changes ^J, . With the appropriate choice of the norms the condition number is the ratio r — of the 
largest to the smallest modulus of the eigenvalues. With this definition, and for the type of system we consider, a 
system can be inverted reasonably if the condition number is smaller than 100 ^ 1000. Note however that a condition 
number can be arbitrarily large but still the system is invertible. This is the case if the RHS of the system is in the 
kernel of the operator. This situation occurs with some preconditioning. 

We now note that, due to the specific form of the Dirac operator Eq. |4j one has 

D{ag, an; cf)^) - I = a {D(g, k; (j)^) - 1) (17) 

where 1 denotes the AN x 47V unity matrix. Since the probability of the (t)^^ does not depend on g and k one is lead 
to introduce the polar coordinates r and 9 of the parameter space (g — rcos{6) n — rsin(0). For a given value of 9 
the spectrum of D evolves straightforwardly : the eigenvectors are then left unchanged and the eigenvalues evolve 
according to 

X''{r,0)^-\>'{ro,e) + l-- (18) 

In a spectral decomposition of D, varying r only changes the relative weights of the eigensubspaces. The value of 
9 fixes the spectrum, and the value of r the relevant part of the spectrum. In general the eigenvalues are complex 
A*^ = + lAf . Let us give a fixed value to 9 and denote the spectrum A'^(r). We choose a reference value rp (one 
can take for example tq = 1) and note A*"' = A'''(ro), one has 

|A'=(r)f = (|aY-2A^, + i) r2 + 2(A^- l)r + l (19) 

So the modulus of the each eigenvalue is a parabola as a function of r. All these parabola intersect at the point 
(r = 0, A = 1). They also intersect each other at others points, and the two extremal eigenvalues change when r 
changes (see Fig. [T]). The eigenvalue labeled by k will reach its smallest value 

for r = |Afc|2^2A'° +1 • Therefore only the eigenvalues with A^ < 1 and Aj' <^ 1 give rise to a small denominator in 
the condition number. When r ~ the eigenvalue of lowest (resp. largest) modulus will be the one with the smallest 
(resp. largest) value of A-^ — 1, therefore the condition number increases continuously from the value 1. In the other 
limit r 1 the eigenvalue of lowest (resp. largest) modulus will be the one with the smallest (resp. largest) value 
of jA'^p — 2Ap^ + 1, and the condition number tends to a finite value (the ratio of the two values above). In the 
intermediate regime, the condition number has a very complicated behavior with a lot of maxima and minima. We 
analyze this behavior is the next subsection for different case. 



This correspond to g = and therefore this is the trivial case of non interacting fermions, it is included for 
illustrating purpose. The evolution of the spectrum of Z) as a function of r = k is straightforward. Performing a 
discrete Fourier transform one finds that the eigenvalues are given by 



Afc = 1-2k^cJ ±2«z /^s2 (21) 
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FIG. 1: Evolution of the square of the modulus of some eigenvalues as a function of r for the free case g = i.e. 6 = 



with Ci, — cos(/ci/) and Si, — sin{ki,). Therefore the condition number behaves as 



c(r — k) — < 
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The condition number diverges at the two values n — ^ corresponding to fc = (0, 0, 0, 0) and n = \ corresponding to 
k = (^,0,0,0) (or permutation when the corresponding are even). The first value | is the critical value, while 
the other is unphysical since it corresponds to a negative mass. 

This is compatible with the framework introduced the introduction of this section. For example Eq [20] becomes 



ruk 



On figure Fig 



zero eigenvalues appear only for r = g and = | 



the evolution of some eigenvalues with r = k is presented. One sees that for 







This case corresponds to k = and describes infinitely heavy fermions. It is unphysical but non trivial. However it 
is instructive to study it from a statistical mechanics point of view, and also because if some continuity is to apply, it 
should not be very different from the 6 small case. In that case the Dirac operator is simply diagonal and the 4x4 
blocs are given by 



(22) 



19 



becomes simply (A'^(r))^ 



Obviously the eigenvalues 1 + g(l)x are all degenerated four times and real. The Eq. 

(1 + r(A'^ — 1))^- Since A'^(r) is linear with r, for any eigenvalue A'^jl there will be a value of r*^ for which A(r*') = 0. 
This is the worst situation since for any fields configuration the determinant of the Dirac operator will exactly vanish 
Y times. This is illustrated on figure 111 where all the eigenvalues as a function of g are shown. The two eigenvalues of 
largest and lowest modulus are emphasized. One clearly sees that the eigenvalue of lowest modulus vanishes for many 
values of r (recall r = g when k = 0). This situation is in contrast with the previous case 6* = f where the eigenvalue 
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of lowest modulus vanishes only twice (for r = ^ and r = j). Here the eigenvalues are "protected" by their imaginary 
part. 

In other words, the (l)^ are N correlated real random variables following the probability distribution Eq. [7j and we 
are evaluating the condition number c(r) which is in that case 

max^^ (23) 
^ ^' mm|l + r(/)j.| ^ ' 

Let us suppose that the have been sorted in ascending order. We note Am the largest negative eigenvalue. Since 
N is very large, we assume m ~ ^ ^^^d there is at least one negative and one positive eigenvalue. The schema on 
figure Fig [2] illustrates the behavior of the spectrum of the Dirac operator for a given (j>x realization. Each eigenvalue 
varies linearly with g. Therefore the condition number is controlled by the eigenvalue of smallest modulus, which is a 
piecewise linear function of g. The selected eigenvalue changes each time g reaches a value gi = — ^ , and reaches 

zero for g* — —-^ for < i < m, which g^ < go < g* < gi ■ ■ ■ ■ Consequently three regimes occur. Firstly when g < g^ 
the condition number is a continuous increasing function of g (homographic) which diverges at go. Secondly in the 
intermediate regime g^ < g < g^ the condition number varies extremely fast diverging m times. Finally for g'!^ < g 
the condition number decreases homographically saturating at a finite value. This is illustrated on Fig. [2] where the 
extreme values gg and g^ are indicated. 

In order to perform analytical evaluation of those tree regimes we simplify the problem by choosing the fields (f>x 
independent with zero mean and a variance given by Eq. [9j It turns out that this simplification does not change 
substantially the average value of the eigenvalue of lowest modulus, as it is illustrated on Fig. |3] This figure shows, 
among other things detailed below, the two curves of the eigenvalues of lowest modulus (curves labeled TV = 131072) as 
a function of r when the 4>x are independent identically distributed Gaussian variables and when they are dependent : 
the two curves are completely indistinguishable. Within this assumption, when the number N of lattice sites increases 
(pm goes to zero as 







Therefore g^ ^ N and the third region shrinks when the lattice size increases. In other words the decreasing of the 
condition number for large values of the coupling constant g at k = is a size effect. On the other limit for small g, 
the first region g < g^ is delimited by the smallest field 0o whose average is given by 



N 



dx (25) 



we see that (^o) diverges extremely slowly with N . To have (0o) of the order of ^, one needs a huge lattice of 

N ^ ^ exp ^ sites. Therefore the first region also disappears in the thermodynamical limit. However this size effect 
will never be seen in an actual computation. Finally we conclude that only the second region survives the large lattice 



volume. Let us recall that in this region and for any fields configuration there are m ~ ^ values of g for which one 
eigenvalue of the Dirac operator is exactly zero. 



In the precedent paragraph the behavior of the condition number for a given configuration of the (p^ has been 
studied. We need now to perform an average over the realization of the (px- For a fixed value of g, different field 
configurations will give very different condition numbers, some of them possibly extremely large. Note however 
that the condition number is not a physical observable, it is only an indicator of how difficult the inversion will be. 
Therefore the most probable value of the condition number is maybe more sensible. From the probability distribution 
of the (px one can easily compute the average of the smallest and largest eigenvalues as a function of g. This is done 
in Appendix |Aj the result is 

(IA.I) - ^/l-xp^ (26) 
(|Aa|) < W21niV (27) 

The eigenvalue of lowest modulus goes to zero as but the prefactor increases extremely fast when g goes to zero. 
Since N goes to infinity first, for any non zero g, goes to zero. The eigenvalue of lowest modulus increases very 
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r 



FIG. 2: Schematic evolution of the modulus of the eigenvalue for one (l>x realization for k — 0. The modulus of the eigenvalue 
of largest and smallest modulus are shown in red. 

slowly with N. This is illustrated in Fig.[3]where we show the eigenvalue of lowest modulus for several lattice volumes. 
On this figure and for iV = 64 32'^ — 131072 we have plotted the result of a "genuine" simulation with dependent 
fields, another simulation with independent fields, a numerical integration of Eq. |A1[ and the approximation Eq. |26| 
The agreement between these four calculation is excellent. We have also plotted the eigenvalue of lowest modulus for 
other values of N to show the size effects. 

The appearance of the three regions described above can be seen in figure |4j On this figures we have plotted the 

average condition number ^ | ^ over 8692 samples as a function of g. If we would have used a smaller discretization 

of g we would have even more sharp peaks. The quantity ^|'^"""'[^ is much smoother since < |Aniin| > never vanishes, 
and also displays the three regimes. Moreover in the quenched approximation it makes sense to consider a particular 
realization since the weight of a consideration does depend only on fi, we therefore have plotted a typical configuration. 
Finally we have also plotted the average condition number without the quenched approximation : this is discussed in 
the next section. 



In conclusion of this subsection, the size effects on this model for k = appear extremely severe: for any fixed value 
of g the occurrence of configurations with arbitrarily small eigenvalues in absolute value grows with the size. This 
is reminiscent of the so-called "exceptional configurations" which have been encountered in the context of quenched 
lattice QCD [12]. 

C. case < f < f 

This region is non trivial since the Dirac operator cannot be diagonalized as in the two previous cases. Nevertheless 
this is where the physics takes place. As it has been done in Ref. [6^, to perform realistic calculation one finds the 
critical line, and one chooses the particular point close to this line where the ratio of the renormalized masses of 
fermion and boson is equal to the physical one. This program has been done successfully giving consistent results for 
g small. However for g around 0.7, the linear system Eq. [TT] becomes ill conditioned, preventing any conclusive result. 

We have computed the condition number for a typical field configuration and the result is presented in Fig. [5j It 
has been shown in IJ] that for small g the critical line, defined as the line where the renormalized mass of the fermion 
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FIG. 3: Smallest modulus of the eigenvalues vs g f or k = and three values of L. For any value of L the numerical integration 
(black) together with the approximation (red) Eq. Al are shown. For the largest size L = 32 x 16"^ the result of a simulation 
averaged over 8692 samples is also shown for both correlated 4>x (magenta) and uncorrelated <^i, (blue). Actually the curves 
are indistinguishable except for small g < 1 when the approximation of the integral is not correct. 



condition number 
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10^ 



100 




quenched average 
quenched typical sample 



unquenched average 



quenched <^a>/ 



0.01 



100 10"^ 10^ 10^ 10^° 



FIG. 4: Condition number vs g for «: = for a 32 x 16^ lattice. The quenched average over 8692 samples is shown in red, and a 
typical quenched sample is shown in blue. For comparison the ratio is also shown in black. The results of the unquenched 
Monte-Carlo simulation is the thick purple line. 
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FIG. 5: Iso-condition number line on a 16 lattice. 



vanishes, is a parabola originating from the point k= = 0. We therefore expected a diverging condition number 
along this line. This is clearly seen on Fig. [5j When the coupling g increases it enters in an ill conditioned region 
where the condition number shows large fluctuations. The localization of this ill conditioned region for k small agrees 
with what we have shown in Sec. |IVB[ Then we can ask if this region where the condition number is small enough is 
not a size effect, as for the case k = 0. 

We see three possible origins for this problematic region. It can be that quenched approximation is not working 
for those values of g. This seems intuitively reasonable since the determinant in Eq |6] precisely give a low weight to 
these configuration with a large condition number. Another possible reason could be the specific choice of the action 
and the discretization of the fermion. Finally there is the possibility that this a fundamental problem of the Yukawa 
model. 



V. THE n = CASE WITHOUT THE QUENCHED APPROXIMATION 



In this section we consider the simple case k = as in the subsection sec. |IVB[ but without the quenched 
approximation. The purpose is to illustrate on this simple case the consequence of the quenched approximation. 
Intuitively the determinant in the probability density Eq. |6] of the (px gives a vanishing weight to the non invertible 
configurations. So we can expect that the configurations to include in the sampling will not have a large condition 
number. But it is possible to have a large determinant, and still a small eigenvalue, for example if one eigenvalue 
is small and all the others are large. These configurations would have a non vanishing weight, but still a very large 
condition number. 

The joint probability of the fields 4>x Eq. [6] can be written as 



n({0j)^exp(-J^ 



(8- 



41n|l 



(28) 



Since this expression cannot be factorized we have written a simple Monte-Carlo algorithm to generate the (j)x^s. We 
use the simple metropolis algorithm [T5]. The normalization factor of Eq. [6]is very difficult to compute, but the ratio 
of the probability of two configurations is very simple to compute (see Eq. [28] ). The Monte-Carlo method use this 
fact to construct a Markov chain which has the desired distribution as a fixed point. In practice, we start from a 
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FIG. 6: <(/!>> in an unquenched simulation on a 32 X 16^ lattice for «: = 



initial (f)x configuration, then we choose at random a site x and try to change the value (px for 4>x + dcf) where d4> 
is a random number normally distributed. We accept this change with the probability min(l, AE') where E is the 
variation of the argument of the exponential in Eq. |28j We do not have a proof that this algorithm converge, the 
difficulty being that the number of states of the Markov chain is infinite. However for all practical purpose it works 
properly if one choose always as a starting distribution for a value g an equilibrium distribution for a close smaller 
value g — Sg. This indicates that the energy landscape is complicated, probably with metastable states. This naive 
algorithm is much simpler than the well known hybrid Monte-Carlo algorithm [T^, but it is sufficient for our purpose 
here. We have compare the two algorithms finding that hybrid Monte-Carlo algorithm is more efficient than the naive 
Monte-Carlo if the parameters are properly chosen, but they both give the same results with a good accuracy. 

Before analyzing the condition number, let us look at the mean value {(fix) (vacuum expectation value). Let us first 
recall that in the quenched approximation, due to the symmetry of Eq. [Tjthe average value {(j)x) is zero. This not the 
case without the unquenched approximation as seen on Fig. |6j even for small g. Indeed performing a ^-expansion of 
Eq. [28] one finds that for any site x 

= 9Y.(^-^yl=o + 0(.9') = 4 + Oid') (29) 
y ^ 

The insert of Fig. shows the slope ^ at the origin : the agreement is very good. Since the average {(j)x) grows with 
g, it seems likely that min(|l -I- g<j)\) will not easily become small. This is indeed confirmed in Fig. [t] where we have 
plotted the eigenvalues of minimum and maximum modulus for both the quenched and unquenched case. It is clearly 
seen that Amin is never small. Finally the average condition number is plotted on Fig. |4] where the drastic effect of 
the quenched approximation is clearly seen : a reduction by six order of magnitude of the condition number. This 
reduction is larger with larger lattice. We conclude that there is no ill conditionned point on this k — line without 
the quenched approximation, whereas it is everywhere ill conditionned in the quenched approximation. 



VI. SUMMARY AND PERSPECTIVES 



We have analyzed in the present paper the appearance of very small eigenvalues of Dirac operator in a Yukawa theory 
with Wilson fermions. The results obtained lead to the conclusion that at finite volume and within the quenched 
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FIG. 7: Amin and Amax in both a quenched and an unquenched simulation for k = (see also Fig. [3| 

approximation, these small eigenvalues are present in an entire region of the phase space. This indicates the existence 
of an ill conditioned region, not just an ill conditioned line, for example the entire k = line is ill conditioned in the 
quenched approximation. Moreover the size effects are exponentially large and consequently a numerical calculation 
can give apparently correct results, which would not survive the infinite volume limit. In other words it does not 
seem possible to determine numerically the ill conditioned region. The origin of this difficulty could be simply the 
choice of the discretization, or it could be the non validity of quenched approximation. This hypothesis is supported 
by a an unquenched calculation for k ~ 0, that is nowhere ill conditioned. But it could also be a problem of the 
Yukawa model itself. Indeed the Yukawa model is not a gauge model and there is no protection against spurious low 
eigenvalues like in QCD [12]. 

In this context we feel that the model should be studied without the quenched approximation. However a boson 
self coupling term has to be added to the Lagrangian to ensure renormalizability. This work is in progress. 



Appendix A: Average of extreme eigenvalues for k = 



In this appendix we show Eq. [26] and Eq. 27 Since the (f>x are normally distributed with zero mean and variance 
given by Eq. |9]the integrated probability distribution of |A| is 

where a is the variance of the (j)x- Then from the definition of the min and after an integration by part, one gets 



(lAminI) 



(l-F|M(x))^dx 



(|A„,ax|) = / I - m^lix))^ dx 

Jo 



(Al) 
(A2) 



introducing ip^y) = 1 — Fc^{y) we have 



(•^min) 



N 
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Since 



then 



(/.(/i) = l-\/^^exp-^/i + 0(/.3) 




</>(|)^^exp(^-.y^^e-=;^^ 

SO 

poo 

NX,nin = / exp 
Jo 

yielding Eq. [24] 

We now study the behavior of Aniax- When N is large the integrand in Eq. |A2| tends to a step function equal to 
one for x < x* and equal to zero for x > x* . One can estimate x* as the unique zero of the second derivative of the 
integrand. Since x* grows when N grows, one can replace x — 1 and x + 1 by a; in the equation '^F\A\{x)^ = 
yielding the equation 

pTT (Xjnax) 1 / (An 

N—\— exp ; ' 



from which Eq. [27] follows. 



2 ga 2 \ ga 



Appendix B: Estimators of boson correlator 



The three estimators give the same correlator 



However the variances are different: 



with 



= {Cl{t)) - {Co{t)f = C\t) + A,{L)C{t) 
- (cut)) - {C\{t)f ^ C\t)+C{0)Cit) 
{cm {C2{t)f = 



2ifet 



3^1, L 



where a is defined in the text Eq. ]9] Consequently one find ctq diverges as L^, tends to a finite value and (t| goes 
to zero as ^ with ~ ^ ^ — t- 0. Only the third estimators C2{t) is self-averaging. 
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